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Abstract. We present a computational approach to the outcoupling in a simple 
one-dimensional atom laser model, the objective being to circumvent mathematical 
difficulties arising from the breakdown of the Born and Markov approximations. The 
approach relies on the discretization of the continuum representing the reservoir of 
output modes, which allows the treatment of arbitrary forms of outcoupling as well as 
the incorporation of non-linear terms in the Hamiltonian, associated with interatomic 
collisions. By considering a single-mode trapped condensate, we study the influence of 
elastic collisions between trapped and free atoms on the quasi steady-state population 
of the trap, as well as the energy distribution and the coherence of the outcoupled 
atoms. 
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1. INTRODUCTION 

The most direct way to create an atom laser, which means an intense directional beam 
of coherent atoms, is based on outcoupling a pre-formed Bose-Einstein Condensate 
(BEC). Key constituents of such a device are the pumping mechanism replenishing the 
trapped condensate as atoms are outcoupled from it and the nature of the outcoupling 
mechanism. This has been demonstrated experimentally(for unpumped condensates) 
leading to the generation of pulsed |JJ El E] and quasi- continuous [H El El HI IHj atom 
lasers. There have been numerous theoretical treatments of the output coupling and 
atom lasers, and a review of the theory of atom lasers can be found in |Hj and references 
therein. Beyond the early rate equation approaches (TUJ HH H2] , these can be crudely 
classified into mean field approaches and treatments based on master equations. Initial 
master equation approaches |IH1 El E3 H3 E| were formulated within the Born- 
Markov approximation, whereas subsequent work [TBI 113 1201 EH 1221 122] has brought 
up issues pertaining to the non-Markovian character of the relevant processes and its 
differences from the optical laser. Moreover, not only the Markov but also the Born 
approximation (standard in quantum optics) is not necessarily valid. Similar issues 
have also arisen in recent years in the context of entirely different phenomena related 
to photonic crystals [21]. In particular, the behavior of an excited atom inside such a 
medium exhibiting a photonic band-gap (PBG) around the atomic transition frequency, 
has been shown to exhibit a number of features also related to the invalidation of 
Born and Markov approximations for the processes involved. As a consequence of the 
unusual density of states appearing in the coupling of the atom to the medium, a 
master equation which is the standard tool in laser theory and quantum optics, can 
not be derived in that case. This rather fundamental limitation has led to the need for 
alternative approaches and techniques, some of which have been proven versatile and 

useful |251 ESI 123 ESI 123 EDI EH- 

Another common approach to the physics of atom lasers is based on mean field 
treatments whereby the dynamics of both trapped and outcoupled components are 
governed by nonlinear Schrodinger equations, with the two systems coupled by an 
external electromagnetic field generating the outcoupling (321 HUH EH ESI ESI EH EH 
I3HJ- Nonetheless, significant insight regarding the breakdown of Born and Markov 
approximations can be obtained by studying the limit of a single-mode condensate, 
as demonstrated in [T%| E3 EH E21 EHj- Although in the limit of no interactions 
between the trapped and the outcoupled atoms, the corresponding master equations 
can be solved exactly, it is still not clear, despite recent work [22J, what is the effect 
of interactions between trapped and outcoupled atoms. In the present work, we show 
that one particular approach [23110111111121, which originated in the context of atomic 
decay in a PBG medium, can be useful in addressing this and other related questions 
pertaining to the outcoupling in atom lasers. The approach we have in mind relies 
on the discretization of continua appearing in the equations of motion that couple the 
few degrees of freedom, usually referred to as "the system", to a reservoir which, by 
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its nature, contains infinitely many degrees of freedom corresponding to a part of the 
Hamiltonian with a continuous spectrum. 

In the following sections, we apply the discretization technique in the context of 
a one- dimensional model for an atom laser, with the limitation to one dimension due 
only to computational constraints. The current paper does not aim at describing the 
physical processes occurring in the recently-observed effectively one-dimensional systems 
113 EH H3 EE3 HZl HEJ which are prone to large phase fluctuations [HU ED] , an issue 
that is dealt with in |51j . Our findings are therefore relevant to and of interest for 
current regimes of three-dimensional atom laser operation [TJI^llIIllllEllBlHIE]- ^ e 
first investigate the validity of the method and study the dynamics of the trapped as 
well as the outcoupled atoms, in the absence of interatomic collisions between trapped 
and outcoupled atoms. We show that, in that context, the method allows us to treat 
practically exactly the problem, even for a range of parameters where it can be treated 
approximately. Then, we extend the treatment beyond the majority of master equation 
approaches that have been addressed so far in the context of atom lasers, by allowing 
interatomic collisions between the trapped and the untrapped (free) atoms. We discuss 
the assumptions under which such processes can be incorporated in our simulations and 
study the system dynamics under their influence. 

2. OUTCOUPLING FREE OF INTERATOMIC COLLISIONS 

We consider a BEC consisting of a large number of bosonic atoms cooled into a single 
eigenmode of a trap. The atoms are coherently coupled out of the trap by the application 
of external electromagnetic fields which induce an atomic transition from the internal 
state (\t)) of the trapped atoms to an untrapped state |/) |S2l E31 ESI EH E2| • Initially, 
only the lowest mode of the trap is populated (condensate mode). We may further 
ignore the population of higher modes even at later times, if the frequency and the 
linewidth of the applied fields are chosen appropriately [THj . 

The dynamics of the system are described by a Hamiltonian of the form 



where H$ and Hr correspond to the trapped (system) and untrapped (reservoir) atoms, 
respectively, while V refers to the output coupling. 

Introducing a set of bosonic operators {a^,a} for the representation of the 
condensate, the corresponding Hamiltonian reads 



where u is the condensate-mode frequency. Accordingly, the part of the Hamiltonian 
for the free atoms is 



H = H S + H R + V } 



(1) 



a 



(2) 




(3) 



where b* k , b^ are the corresponding bosonic operators, with k being the momentum of 
the center-of-mass motion. 
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2.1. Output coupling 

The condensate is coupled to the reservoir through the interaction term 

V = h f + °° dkg{k){b\a + b k a)). (4) 



In general, the form of the coupling constant g(k) is determined by the type of 
mechanism applied to coherently couple the trapped atoms out of the trap (output 
coupler) [3H]. For a Raman output coupling mechanism, the outcoupled state mimics 
the trapped condensate at short times (Hill EH]. Assuming Gaussian like profile for the 
ground state of the harmonic trapJTHJ [TH], the coupling constant can be approximated 
by 

9(k) = j^yT^M-k 2 /^ (5) 

where T and a denote the strength and the width of the coupling, respectively. 
Additionally, the frequency u>k is related to the corresponding momentum by Uk = 
hk 2 /2m. It is this quadratic dependence of u>k on k, which together with the functional 
dependence of g(k) on k, give rise to mathematical difficulties even in this simple model. 
Note that, the dependence of uj^ on k is in fact identical to the analogous dispersion 
relation in PBG materials, which is responsible for the mathematical difficulties in that 
context. 

For the one-dimensional problem under consideration, the density of states which 
are available to a free atom can be determined by the dispersion relation as follows 

dk I m 

e(w), (6) 



duo 



2hu 

where m is the atomic mass and the usual step function. Taking advantage of the 
symmetrical shape of the coupling and the even parity of uj^, we may reduce the fc-space 
only to the k > sub-space. The spectral response of the output coupling is then of 
the form 

V G (u) = 2|^(^)| 2 pM = [ ' J 6(o;), (7) 

7T \ UJ 



where the effective coupling constant C is given by 

C = r</^ and a = — . (8) 
V a m 

As a and T tend to infinity, we obtain the broad-band limit of equation (|7|). 

v bb (oj) = -^en, (9) 

7T \/0J 



which is identical to the corresponding spectral response of PBG continua. Note that 
our model does not consider any losses due to gravity or Raman-momentum kick which 
could, however, have been included. 

Given the Hamiltonian (0) one may proceed to derive the Heisenberg equations of 
motion for the operators of interest, examples of which are discussed below. At this 
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point, let us recall that in problems where the Born and Markov approximations are 
valid, as for example an excited atom decaying spontaneously in the vacuum of open 
space, the reservoir (continuum) variables can be eliminated. This procedure leads to a 
Markovian master equation governing the evolution of the system. This is what can not 
always be done in the present context. In the following subsection we present, for the 
model Hamiltonian under consideration, an approach capable of providing not only the 
time dependence of the number of atoms in the condensate, but also the distribution of 
the outcoupled atoms in frequency domain, irrespective of the strength of the coupling 
between the system and the reservoir. Moreover, it is applicable to more general coupling 
shapes and suitable for the inclusion of additional interaction terms in the Hamiltonian 
(see section EJ) for a more realistic description of the physical situation. 



2.2. Discretization of the reservoir 

In order to deal with the structured reservoir, we follow the discretization approach 
developed in the context of PBG continua [23 HOI EH Q2j Briefly, we substitute the 
reservoir for frequencies within a range around ujq {uji ow < uj < u up ) , by a number (N) of 
discrete modes, while the rest of the atom-mode density is treated perturbatively since 
it is far from resonance. In general, there is no unique way of discretizing a continuum. 
For the sake of illustration, in the present work we apply two different discretization 
schemes for the Gaussian-like profile and its broad-band limit, respectively. 

Specifically, for the Gaussian-like coupling we adopt a uniform discretization 
scheme, choosing the frequencies of the modes to be 

Wj = ui ow + jSco, (10) 

where the mode spacing Sou, is determined by the upper-limit condition of the 
discretization, namely 

u uv = ui ow + N5uj. (11) 

The corresponding coupling for the j mode, is determined by the spectral response ((Zj) 
as follows 

Q) = V G (uj)8w. (12) 

Alternatively, for the broad-band coupling (equation EJ), we may choose the integral 
form of equation (fT2"j) . obtaining thus frequency-independent couplings for the modes, 
determined by 

G 2 N= / duV BB {u), (13) 

while the frequency for the j reservoir atom-mode is given by 

= u} iow + j 2 5u, (14) 
with the spacing Sou determined by 

uJ up = Wow + N 2 5uj. (15) 
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Finally, note that in A;-space the discretization would rely on the following equations 

kj = h ow + jdk, Qj = g{kj)5k, k up = k low + j5k, (16) 

from which we may obtain equations (jlOjl - lfTo'j) . using the definition of the density of 
atomic states (equation EJ). 

2. 3. Heisenberg equations of motion 

In the Heisenberg picture, the evolution of an arbitrary operator A is governed by 

£ = 4^ (17) 

where H is the Hamiltonian, which for the problem at hand is determined by equations 
©-©• Thus, for the operators of the system and the reservoir we obtain the following 
equations of motion 
da 

— = — iui^a — 2% / dkg(k)b k dk, (18) 
at Jo 

— = - iu k b k - ig(k)a, (19) 

= - 2% (°° dkg(k)(a% - bta). (20) 
dt Jo 

The solution of the above differential equations can be obtained in terms of 
the inverse Laplace transform. Furthermore, as has been noted in [T%1 IT9"1 |2*U] . the 
solution for the mean number of atoms in the condensate ((a'(t)a(i)}) can be written 
as (a*(t)a(t)) = (a* (t)) (a(t)) , where (a(t)) can be obtained as solution of equation (JTHJ), 
assuming; the initial condition (a(0)) = ^/(at(O)a(O)), while (a f (t)) = (a(t))*. This 
assumption stems from the fact that the initial state of our model corresponds to a 
BEC in the atomic trap, which is a coherent state \f3) of definite global phase and 
thus (^(0)0(0)) = |/3(0)| 2 = (a(0))*(a(0)). It is the bilinear form of the Hamiltonian, 
however, that preserves this coherence in time. What has to be noted here is that, 
irrespective of the initial conditions, we may obtain the evolution of {a){t)a{t)) by 
combining equation ()20|) with the following equations of motion for a\t)bk{t) and 
b\(t)b q (t) respectively: 



d(a ] b k 

dt 
d(b\b q 



poo 

i(u k - u )a^b k - ig k a ] a + 2i dqg(q)b\b k , (21) 
u k )b\b q + ig(k)a ] b q - ig{q)ab\. (22) 



dt 

Discretizing the continuum over a range of frequencies (ui ow < uj < u up ), the 
integrals are converted to sums over discrete modes of frequencies 0Jj{ m )i with j(m) 
running from 1 to N, while the number of discrete modes covers only part of the 
reservoir. Assuming, as must be the case, that the discretized range is sufficiently 
large to incorporate the essential effects of the system-reservoir coupling, the remaining 
part of the continuum (uj < ui ow and u > oj up ) can be eliminated adiabatically to 
order of \g(k)\ 2 , in the standard fashion. Specifically, the derivative of b k in equation 
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(Unj) for Uk outside the discretized range, is set equal to zero, yielding the approximate 
solution 

h(t) ~ a(t), for to k < u low and to k > uj up . (23) 

to k 

These approximate solutions are then fed back into the equations for the operators in 
the discretized part of the continuum. As a consequence of this adiabatic elimination, a 
shift appears in the equations involving derivatives of the operators a(at) and a^b k (ab\). 
Furthermore, since the equations of motion (|2U|) - (|22|) constitute a closed set of 
differential equations we may deal with expectation values of the operators belonging 
to the discretized space of the continuum, obtaining 

,/<(/ '' f) =2j:G^{(a%)}, (24) 



dt 

d(a^bj) 

dt 
d(blbj) 



N 

- i(uj -uj + S)(a 1 b j ) - iQ j (a ( a) + i ^ Qm^lrfij), (25) 

m=l 

- i(ujj - Urnj^bj) + iQ m {a)bj) - i<3j(a ] b m )*, (26) 



dt 

where for simplicity of notation, u r , Q T and b r stand for to kr , Qk r and b kr respectively, 
with r e {j,m}, while denotes the imaginary part of the expression inside the 

curly braces. The shift term is given by 

S= / du^-+ / duo^-, (27) 

JO CO Ju>up tO 

where V(uj) = T>g(bb) and the lower bound of the first integral has been set to zero as 
dictated by the step function in the spectral responses of the reservoir. Furthermore, in 
the above set of equations we have used the fact that for each pair of arbitrary operators 
A and B we have (AB) ] = B^A ] and thus (AB)* = The same equality allows 

us to deal with half of the matrix elements ifr m bj), those which involve m < j since the 
elements involving j < m can be obtained by (b^bj) = {b^b m )* . 

In the following sections, we present some selected results for the above discretized 
set of differential equations. Specifically, for various values of the parameters and for 
both Gaussian and broad-band couplings, we investigate the time dependence of the 
number of atoms in the condensate and the energy distribution of the outcoupled atoms 
(blh). 

At this point, some general remarks need to be made concerning the expected 
accuracy of the discretization approach. Two main parameters enter the system of 
equations, namely the number of modes N and the upper(lower) limit of the discretized 
part of the continuum. The number of modes will determine the maximum time for 
which the numerical solution of the discretized system can be considered accurate. This 
is already known from the experience gained in PBG media. If one wants to obtain a 
correct solution for longer times (and this is the case for larger values of coupling, when 
a steady state is reached more slowly), one has to increase the number of modes and 
the number of differential equations, which scales roughly as iV 2 . 
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On the other hand, the parameters u up and ui ow must be chosen so that u up — lui ow 
is large enough to include all modes with significant role (for sure it must be chosen such 
that ujiow < u < u up ) but also such that the spacing 5u will produce a smooth enough 
energy spectrum. As a general remark pertaining to the results discussed below, we 
note that the convergence of our calculations has been checked in terms of the number 
of discrete modes as well as the range of discretization. As is depicted in figure [TJ 
in order to reproduce the results obtained by means of the Laplace transform in the 
context of Gaussian couplings(solid line) [20J, we find that we need at least 500 discrete 
modes (dot-dashed line) and a value of 3 for u up in units of C 2 / 3 . As far as the short- 
time and the long-time behavior are concerned (in which we are mainly interested), we 
may obtain the dynamics of the system even for 50 or 60 modes and u up = 3 (dotted 
line), avoiding thus time-consuming numerical calculations which are expected to add 
quantitative and not qualitative corrections to our results. For the broad-band coupling 
(inset) the situation is slightly different, in the sense that the solution obtained with 50 
discrete modes (dotted line) is almost indistinguishable from the exact solution (solid 
line). Nevertheless, as a consequence of the discretized continuum, it exhibits revivals 
for longer times, which are expected to appear for later and later times as we increase 
the number of discrete modes. 

To facilitate comparison with the results of Moy et al. 1201, throughout this 
paper we consider m = 5 x 10~ 26 kg, ujq = 2n x 123SCC" 1 and a = 10 6 m _1 . Since, 
however, the time scale on which the main dynamics take place is expected to be 
strongly dependent on the coupling parameters, from now on we present the dynamics 
of the system as a function of the dimensionless time C 2 ^t. 

2-4- Dynamics 

In figure |2fa), we present the time evolution of the normalized mean number of atoms 
(population) in the condensate mode for Gaussian coupling and various values of T. As 
was expected, for weak couplings in relation to u , the trap mode exhibits a Markovian 
behavior in the sense that its population decays exponentially to zero (long-dashed line). 
As we increase the coupling strength, however, the trap population begins exhibiting 
non-Markovian features. Specifically, after a transient regime where part of the initial 
population is lost, the system undergoes Rabi-like oscillations. These oscillations can 
be interpreted as a beating of the system between two different paths associated with 
the two characteristic frequencies of the problem, namely the trap-mode frequency 
cuo and the singularity oj = 0, where the spectral response is peaked. The crucial 
role of these two characteristic points in the system dynamics, is also reflected in the 
energy distribution of the outcoupled atoms [figures GJb) -(d)] which, for relatively large 
couplings, exhibits peaks at uj = u and u> ~ 0, respectively. In the weak-coupling 
regime, the peak at u>o is prominent. As we enter the strong coupling regime, however, 
a larger part of the continuum is involved and thus, on the one hand the distribution 
becomes broader while on the other, the Gaussian coupling tends to couple the cavity 
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mode not only to modes with frequencies in the vicinity of uj , but also to low-energy 
modes (w « 0). Furthermore, the central frequency of the distribution is shifted towards 
higher frequencies. This behavior is analogous to the shift in the levels of a two-level 
atom coupled to a photonic continuum and in the context of atom lasers has also been 
noted by Jeffers and co-workers 

From the physical point of view, the oscillations in the trap population stem from 
the outcoupled atoms which are fed back into the trapped state, reexciting thus the 
condensate mode. Due to this back action of the reservoir, the system finally reaches 
a steady state where the condensate mode is partly excited. With increasing coupling 
strength, the oscillations become more pronounced and faster, while the population 
trapping in the cavity increases. As is evident in figure El the non-zero steady-state 
trap population is also associated with strong coherence between the atoms in the trap 
and the outcoupled atoms (upper graphs), as well as between the free atoms (lower 
graphs). Although for weak couplings the coherence is negligible and is restricted to 
the vicinity of the trap-mode frequency, it increases rapidly for increasing values of the 
coupling strength, while simultaneously it spreads out to a larger number of modes. This 
behavior clearly indicates the departure of the system from the Markovian dynamics. 
In figure El we note that the index of the discrete mode corresponding to the trap- mode 
frequency u , is given by 



where [x] denotes the integer part of x. 

For the broad-band coupling (figure HJ), the dynamics are basically the same. 
The broad-band nature of the coupling, however, allows a significant interaction with 
more continuum states. Thus, the system reaches steady state much faster, while the 
oscillations in the cavity population are not so pronounced. Moreover, the energy 
distribution of the outcoupled atoms is much broader than in the Gaussian coupling. 

3. INTERATOMIC COLLISIONS 

Up to now we have not considered any collisional interactions between trapped and free 
atoms in our model. As a result, we had to deal with a Hamiltonian bilinear in the field 
operators and thus with a closed set of linear differential equations for the operators of 
interest. In real condensates, however, interatomic collisions take place and affect, in 
most cases, the dynamics of the system significantly. Specifically, considering a small 
rate of output from the trap, the output atomic beam is dilute and thus, we may neglect 
collisional interactions between free atoms. We may further ignore inelastic collisions 
between untrapped and trapped atoms, since the corresponding mean free path is much 
larger than the dimensions of a typical atomic condensate in the trap. 

In our subsequent discussion, we focus on the elastic collisions of the untrapped 
atoms with those in the trap. Considering low-energy s-wave scattering of length 



(uj - u low )N 




Jo = 



^up ^low 
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a s between trapped and free atoms, we may adopt a standard hard-core interatomic 
potential which for a three-dimensional model in position representation is 

U tf (r-r') = U$ D) d(r-r'), (29) 

where U^ D ^ = Arc?i 2 a s /m. For the model at hand, however, we have to derive an 
one-dimensional, approximate form of equation (J29|) . capable of providing us with an 
estimate of the importance of interactions in our model. In the case of tight transverse 
confinement, the scattering strength can be well approximated by integrating out the 
transverse directions leading to [HI] 

U tf (x - x') « - 1L ^5{x - x 1 ) = 2ha s to ± 5(x - x'), (30) 
2 ™ ho 



where u± the transverse confining frequency and a^ = ^h/mu± the corresponding 
harmonic-oscillator length. In this paper we are not concerned with the limit of reduced 
dimensionality and we have only limited ourselves to one- dimensional simulations 
for computational ease. The corresponding part of the Hamiltonian then takes the form 

/+oo r+oo 
dk / dqn(k,q)b k b q a'a, (31) 
-oo J — oo 

where, for a Gaussian ground state of the trap, the coupling constant n(k, q) reads, 



k(/c, q) = Na s u± exp 



(k - qf 



(32) 



8a 2 

with M being the number of trapped atoms. Note that the coupling constant is 
symmetric with respect k and q and thus Vtj is Hermitian. 

In analogy to the output coupling, the spectral response for interparticle 
interactions can be defined as 

M G (u,u') = 4\ K (u,u')\ 2 p(u)p(u') = C M eXp[ ~ (v ^ / I^ 7)2/2a] eM9(^), (33) 
where the effective interatomic coupling Cm is 

C U = (34, 
As a and T tend to infinity, we obtain the broad-band limit of equation ()33)1 : 

Mbb(co,uj') = C M -jL=Q(u)e(u'). (35) 

For a s ~ lnm and m ~ 5 x 10~ 26 kg, we have Cm ~ I0~ 9 (j\fu)±) 2 sec. A reasonable 
estimate of realistic interatomic coupling strengths in our model, can be obtained by 
noting that, for Af ~ 10 3 atoms and u± ~ 2n x lOOOsec -1 , a typical value for Cm is 
C M ~ lO^ec" 1 , or else C M ~ IOC 2 / 3 . 
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3.1. Heisenberg equations of motion 

Having now included interatomic collisions in our model, the Hamiltonian is not bilinear 
and the equations of motion for the operators of interest are no longer linear. Specifically, 
using equation (|T7|) we obtain the following set of differential equations 

dt 
d(a^b k ) 



2i / dkg(k)(a f b k -bla), (36) 



o 



-i[uj k - Uo)a h 



/•OO /"OO 

ig k a}a + 2i / dqg(q)b\b k — 2i \ dqK(k ) q)a)b q a) c 
Jo Jo 



dt Jo 

/■oo roo 

+ 4z/ dq dq'K(q,q')blb q >a%, (37) 



^ = —i(uj q — uj k )b\b q + ig(k)a)b q — ig{q)ab\ — 2i f dk'n(q, k')a)ab\b k i 

Ctv J 

roo 

+ 2i / dk'K(k\k)a)ab\,b q , (38) 



o 

which can not be solved by means of the Laplace transform. 

In order to proceed to the solution of the above equations in the context of the 
discretization approach, we have to incorporate appropriately the interatomic collisions 
in the discretization schemes developed earlier in sectional To this end, we assume that 
the interatomic coupling spans the same discretized space with the output coupling and 
thus does not affect the far off-resonant b k operators. This assumption, on the one hand, 
allows us to adiabatically eliminate the same b k as before, while on the other to use the 
same discrete modes for both the output and the interatomic couplings. Note that in 
the case of trap-trap interatomic interactions, an adiabatic elimination procedure of the 
high-lying modes has been shown to lead to the renormalization of the (bare) interatomic 
potential to an effective one given by the usual pseudopotential of equation (29) [57)1 loT)] . 
We believe it is therefore consistent here to adiabatically eliminate high-lying modes 
(of the output spectrum) to obtain s-wave scattering between trapped and outcoupled 
atoms. Consideration of low-energy s-wave scattering for interatomic collisions implies 
that our discretization should be applied under the constrain k up a s << 1. To confirm 
this, note that for the parameters used in the paper, k up w 10 6 m _1 and a s w 10~ 9 m and 
thus k up a s ps 10~ 3 << I. 

Extending equations (fT^j) and (|T3*|) to the two-dimensional (k m , kf) space we 
introduce the following coupling constants which involve pairs of modes: 

GL =M G {uj m ,u j )8uj5uj. (39) 



J2 A r2 
*int 



- 1/ p 



N 2 = duo' dujM BB {uj,uj'), (40) 



and correspond to the uniform and the non-uniform discretization scheme, respectively. 
We thus obtain the following set of equations of motion for the expectation values of 
the operators involving the discrete part of the continuum 

^ = 2E£,3{(at6,)}, (41) 
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dialb) N N 

— —2- = -i(ujj -uj + S)(aJbj) - ig j (a f a) + i E GmiKfij) - i E GjmWb m a^a) 

at m=l m=l 



N 

E 

Lm=l 



+ i E Gim(bjb m a%), (42) 



fj/hi h.\ N 

Aj^lL = _Hu. _ cuJiblb,) + ig m (a%) - iQMbmY -ijy^ablbt) 

N 

+ iT,Gi m (a j ab]b j ), (43) 
i=i 

which is not closed. The structure of these equations is a typical case of equations of 
motion emerging from Hamiltonians involving terms of order higher than bilinear. The 
main feature is the presence in the right hand side of operators of order higher than 
the one whose derivative is considered. Consideration of differential equations for those 
higher order correlation functions leads to the appearance of terms of even higher order 
and so on. The system of equations, in other words, does not close and there are no 
general exact remedies. 

One way to obtain an approximate solution is to decorrelate the higher order 
correlation functions into products of lower ones. Again, there usually is more than 
one way of doing so, the only guide being the demand for a closed system of equations. 
Perhaps the most straightforward decorrelation suggesting itself in the system of 
equations (|4*T]) - (f4*H|) is to decorrelate products of four operators in a way such that 
the decorrelated parts involve equal number of raising and lowering operators i.e., 
{a^bmO^a) « (a t a)(a t 6 m ), (b\b m a)bj) (b\b m )(a)bj) and (a^ab^bi) « (a 1 a) (b^bi) . We 
thus have the set of equations A: 

^ = 2E^3{(^}}, (44) 
d(a)b) N N 

at m=l m=l 

N 

+ i(a%) E Sim(b}b m ), (45) 

Z, 771=1 

dltf h ) N 
ai l=l 

N 

+ i(<Ja)Y J Gim(b\b 3 ). (46) 
i=i 

Applying the commutation relation [b m , b\] = 5 mn in equation (J4l?|) and decorrelating in 
the same fashion the resulting equation, we have an alternative form for equation (j4T)j) : 

MM = _ a; w )(6j n 6 i ) + iQ m (a\) - iQjl^KY ~ i(a j b m )* E Gfl{a\) 

at l=l 

N 

+ i(a%)Y,Gim(a%y, (47) 
i=i 
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which combined with equations (|44j) and (|45|) constitute another closed set of differential 
equations (set B). 

3.2. Dynamics 

In the absence of interatomic collisions the effective coupling constant C is the one that 
determines the time scale of the initial transient regime i.e., the speed at which the 
trap-mode loses population for short time. In the present case, however, interatomic 
collisions are expected to be the dominant processes for short times, determining thus the 
corresponding time scale. Both sets of equations (A and B) predict the same dynamics 
for the system. In figure we plot the solutions for the population in the trap as 
a function of the dimensionless time C 2 ^t obtained through the set of equations B. 
Clearly, the rate at which the trap loses population decreases with increasing values of 
Cm while the trap-mode excitation undergoes oscillations which reflect the interference 
between the two possible decay routes for the system, corresponding to ojq and u> = 0, 
respectively. 

The functional dependence of JA(u, to') on uj and u', respectively, reveals the main 
role of interatomic collisions, namely the decoherence. We know already (section |2J) 
that the coherence between the atoms is associated with the non-Markovian behavior of 
the system and the partial depletion of the excitation of the initially highly populated 
condensate-mode. Under the influence of interatomic collisions, however, the coherence 
between the outcoupled atoms is reduced and is restricted to the vicinity of the trap- 
mode frequency (see lower graphs in figure leading thus to a significant decay of the 
initially trapped atomic population, as well as a narrow distribution of the energies of 
the outcoupled atoms around ojq [figures Efa)-(d)]. This is in analogy to the behavior 
of the system in the absence of interatomic collisions and for weak output couplings 
(see lower graphs in figure El for Y = 10 4 — 10 5 sec -2 ). It could thus be argued that the 
dynamics of the system are governed by a set of equations of motion of the form 

^ = 2j2gM(*%)l (48) 

at j=1 

d ^ a ^ = -i(ujj - uj + S)(a^bj) - iQ j {a [ a) + iQj{b]bf) - i{o)a) ^ Qjm(a^b m ) 
dt m=l 

N 

+ i{a%)J2Smm{bl n b m ), (49) 

771=1 

rMtf h \ N 

= -2S m 3{(at& m )} - 2j2QirrM(^b m )(a%y}, (50) 

where no coherences between the outcoupled atoms are involved. As depicted in the inset 
of figure El this argument is indeed true for short times and relatively large interatomic 
couplings, where interatomic collisions dominate over the output coupling. As was 
expected, however, the above set of equations fail to describe the decoherence in our 
model for larger times. Alternative decorrelation schemes and their predictions about 
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4. SUMMARY 

We have presented an alternative approach to a simple model for atom laser outcoupling 
which is capable of circumventing certain mathematical difficulties arising from the 
breakdown of the Born and Markov approximations, inherent in essentially any model 
of the process. The approach is computational, relying on the discretization of the 
continuum representing the reservoir of the outcoupled atoms and has been recently 
developed in the context of PBG continua. To provide a calibration of the possibilities 
offered by the method, we presented representative results, some of which serve as 
a basis for comparison with previous work relying on different methods, while others 
demonstrate how this approach can go beyond the above mentioned limitations. 

In particular, we have studied the dynamics of a single-mode trapped condensate 
which, in order to form an atom laser, is outcoupled by both Gaussian and broad- 
band couplings. In the model considered here, there is no mechanism to remove the 
outcoupled atoms, and the ensuing coupling between trapped and outcoupled atoms 
provides (in certain limits) a striking manifestation of the non-markovian nature of the 
outcoupling. For example, considering initially a situation in which condensate and 
outcoupled atoms do not interact, we have shown that, after an initial transient regime 
where part of the trapped atomic population is lost, the system reaches a steady state for 
which the population inside the trap has not been totally depleted. This non-markovian 
nature is most pronounced when outcoupling at a rate larger than (or comparable to) 
the trap frequency, and is also evident in the coherences between trapped-outcoupled, 
or outcoupled-outcoupled atoms, which additionally experience a peak around uj ~ 0. 

These coherences thus play a significant role in the time evolution of the system. 
Perhaps more significant is the effect that the addition of interactions between 
condensate and outcoupled atoms has on the outcoupled spectrum. It has been 
previously argued that interatomic interactions could potentially destroy this 'bound 
mode' of atoms within the trap. In this work we have explicitly considered the effect 
of interactions between trapped and outcoupled atoms and have shown directly that 
realistic interaction strengths will largely destroy this 'bound mode'. Interatomic 
interactions also have a drastic effects on the coherences between trapped and outcoupled 
atoms. We have thus shown that our computational approach is capable of exploring the 
entire domain of outcoupling from non-markovian to markovian, and can additionally 
include atom-atom interactions. 

It should be stressed that this paper only describes the mechanism of outcoupling, 
and does not deal with the full complexity of an actual atom laser. In particular, 
we have assumed that we begin with a single-mode condensate, without any thermal 
excitations being present. Thus, strictly-speaking, our treatment is restricted to the 
zero temperature limit. Consideration of excitations into our model will lead firstly, 
to a decrease in the amount of coherence present in the trapped system. This will 
subsequently limit the coherence transferred to the outcoupled beam (atom laser), with 
the extent of the decrease from the purely condensed case, depending on the particular 
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choice of the outcoupling, since some thermal component will be inevitably outcoupled. 
This point has been discussed in |33J ES], whereas a three-dimensional treatment for 
Raman pulsed atom lasers has been recently performed [H7j . 

The model considered here, does not account for a mechanism of removal of the 
outcoupled atoms. Such a mechanism arises naturally from the effect of gravity, and 
can also be controlled by the amount of momentum imparted to the outcoupled atoms 
via Raman outcoupling. When the outcoupled atoms are removed from the system, 
the only way to reach the desired steady state is by continuously pumping the trapped 
condensate. Both of these features are essential in the steady-state operation of an atom 
laser, and can be included without particular difficulty in an extended version of our 
presented model. However, it bears repeating that the main aim of this paper was not 
to describe the full atom laser operation, but rather to introduce an alternative method 
for describing the process of outcoupling. A model dealing with a cw atom laser in a 
one-dimensional waveguide by means of a Langevin equation, where issues such as the 
laser linewidth can be addressed has been discussed by one of us elsewhere |HTj . 

One of the important features of the method employed in this paper is the versatility 
in handling any form of functional dependence of the outcoupling on the atomic 
momentum. Furthermore, it is amenable to the study of further properties of the 
system such as higher order correlation functions, albeit, as it appears for the moment, 
within the limits of some form of decorrelation. Although the various decorrelation 
approximations we had to resort to have given more or less qualitatively compatible 
results, still the task of narrowing the possible choices needs to be addressed. One avenue 
in that direction would be the examination of higher order of correlation functions, with 
the ultimate limitation being the size of the computation, as the model becomes more 
complex. Computational possibilities, on the other hand, are not frozen but increase 
constantly. Clearly, analytically solvable models are indispensable in the insight they 
offer, but computational approaches are often necessary in examining details of the 
dynamics of a complex system. 
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Appendix A. Alternative decorrelation schemes 

In general there is no unique way of decorrelating a system of coupled Heisenberg 
equations of motion. In the main body of the paper, we have already presented two 
different decorrelation schemes. We discuss two more in this appendix. Applying the 
commutation relation [6^,6^-] = S m j in equation (JI2J), we obtain 

d(a>b) N N 

— — — = -i(cOj -uj + S)(a%) - iQj(a^a) + i ^ Gm{vU>j) - % E Gjm(a)b m a ] a) 

m =l m= l 

N N 

+ i Gimibmbja^bj) - i{a%)^Gii (A.l) 

l,m=l 1=1 

which can be decorrelated as follows 

dlatb-) N N 

— — = -i(ujj -uj + S)(a' i b j ) - iGj(a^a) + i GmiKfij) ~ ^( flta ) E Gj m (a^b m ) 

m=l m=l 

N N 

+ t E Gi m {cJb m ){b\b 3 ) -i^b^Gu- (A.2) 

l,m=l 1=1 

This equation of motion can replace equation (|35Jl in the sets of equations we presented 
in the previous sections. The resulting sets of equations, however, predict substantially 
different behavior for the system than the one discussed up to now. Specifically, we have 
found that the number of the outcoupled atoms decreases significantly for large values 
of Cm in relation to C 2 / 3 . Thus, the atomic population of the trap tends to remain 
trapped forever ((a)(t)a(t)) ~ 1). Actually, this behavior stems from the last factor 
in equation (jA.2j) . which leads to a non-dissipative Rabi-like oscillation of the cavity 
population close to unity i.e., {a\t)a(t)) ~ 1 — exp(— iJ2i Gut). We may overcome this 
problem by decorrelating the last term in equation as (b[b m a%j) ~ (a)b m ) (bjbj) 
obtaining thus equation (|A.2|) without the term involving Gu- 

Using the corresponding commutation relation for the cavity-mode operators i.e., 
[a, a 1 ] = 1, we obtain from equation (jlHjl 

d(tf b) n n 

' "' " =^(u j -u m )(blb j )+ig m (a%)-ig j (a^b m )*^ 



dt 

N N 



1=1 



~ i^Gjiiaa^blfii) + i^Gimiaatbjbj), (A. 3) 

i=i i=i 

which after decorrelation reads 

dltf b) n n 

AjrJL = _i( u . _ uJiblbj) + iG m (a%) - iGj(a j b m )* + iJ^Gji (blk) - i^imibjbj) 
dt i=i i=i 

N N 

~ iY.$iMb m Y{a\) +i^Gi m {a%)*{a%). (A.4) 
i=i i=i 

Although this equation of motion, accompanied by equations (jSj), (JEJ), predict mainly 

the same short-time dynamics that we have discussed in sectional it overestimates the 

quasi-steady state of the trap-mode as well as the oscillations in the cavity excitation. 
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In other words, equation (jA.4J) seems to overestimate the coherences in our system, due 
to the appearance of terms like J2iLi Gim(bjbj) in its right-hand side. 
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Figure captions 




Figure 1. The normalized mean number of atoms in the condensate as a function 
of the dimensionless time C 2 ^ 3 t for Gaussian (r = 10 6 sec~ 2 ) and broad-band (inset) 
couplings. The solid lines correspond to the exact solutions obtained by means of the 
Laplace transform. We also plot the discretization solutions for N — 50, uoi ow — and 
lu up = 3.0C 2 / 3 (dotted line) as well as for N = 1000, u low = and uj up = 3.0C 2 / 3 
(dot-dashed line), respectively. Inset: N — 50, uoi ow = and uj up = 10. 0C 2 ^ 3 (dotted 
line) 
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Figure 2. The evolution of the normalized trap population, as a function of the 
dimensionless time C 2 ' 3 t (a) and the energy distribution of the outcoupled atoms at 
t = 30C~ 2 / 3 [(b)-(d)], for Gaussian coupling and various coupling strengths. The 
calculations are for 1000 discrete modes, lui ow — and u up — 18C 2//3 (r = 10 4 sec~ 2 ), 
uo up = 6C 2 / 3 (r = 10 5 sec- 2 ), lo up = 3C 2 / 3 (T = 5 x 10 5 - 5 x 10 6 se C - 2 ). 
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Figure 3. The magnitude of the matrix elements (a' (t)bj(t)} (upper) and (b' m {t)bj{t)) 
(lower) for Gaussian coupling at C 2 ' 3 t = 30 and various coupling strengths T. For the 
lower graphs, dark and white regions represent matrix elements with negligible and 
large magnitude respectively. The calculations are for 60 discrete modes and cutoff 
frequencies: u low = 10C 2 / 3 ,lu up = 18C 2 / 3 (r = 10 4 se C - 2 ); u low = 0,cj up = 6C 2 / 3 
(r = 10 5 sec- 2 ); uj iow = 0,u up = 3C 2 / 3 (f = 5 x 10 5 - 10 6 se C - 2 ). 




Figure 4. The evolution of the nomalized trap population, as a function of the 
dimensionless time C 2 ^ 3 t (upper graph) and the energy distribution of the outcoupled 
atoms at t = 15C~ 2 / 3 (lower graph), for broad-band coupling and various coupling 
strengths. The calculations are for 100 discrete modes, wi ow — and uj up = 20C 2 / 3 . 




Figure 5. Effect of trap-free interatomic interactions on the time-evolution of 
the normalized trap population for Gaussian (a) and broad-band (b) couplings. All 
calculations have been obtained in the context of set B. (a) The coupling constant is 
r = 10 6 sec~ 2 and the population is plotted for various interatomic-coupling strengths: 
C M = (long-dashed line), C M = 4C 2 / 3 (solid line), C M = 8C 2 / 3 (dotted line) 
and Cm — 12C 2 / 3 (dot-dashed line). Discretization parameters: N — 60, u>i ow = 0, 
u> up — 3C 2 / 3 . Inset: The solid and the dotted lines correspond to a propagation of set 
B and equations (gEJ>-|HiJ) respectively, for T = lOSec^ 2 and C M = 12C 2 ^ 3 . (b) The 
solid line is for Cm — 0, the dashed line for Cm — 1C 2 / 3 and the dot-dashed line for 
C M = 2C 2/3 . Discretization parameters: N = 60, uji ow = 0, 5ui = 0.004C 2/3 . 
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Figure 6. Effect of trap-free interatomic interactions on the energy distribution 
of the outcoupled atoms (upper graphs) and the magnitude of the matrix elements 
(b} n (t)bj(t)) (lower graphs). All the graphs are for Gaussian coupling (r = 10 6 sec~ 2 ) 
and correspond to C 2 ^ 3 t = 25 and various interatomic-coupling strengths. Dark and 
white regions in the lower graphs represent matrix elements with negligible and large 
magnitude respectively. Discretization parameters as in figure 



